In [1]:
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[1]:
In [2]:
import matplotlib.pyplot as plt
import numpy as np
import healpy as hp #1.40
import matplotlib.style
import matplotlib as mpl
mpl.style.use('classic')
from pylab import arange, show, cm
cmap = cm.bwr
cmap.set_under('w')

%matplotlib inline
params = {
	'axes.labelsize': 24,
	'font.size': 24,
	'legend.fontsize': 20,
	'xtick.labelsize': 22,
	'ytick.labelsize': 22,
	'lines.linewidth': 4,
	'text.usetex': False,
	# 'figure.autolayout': True,
	'ytick.right': True,
	'xtick.top': True,

	'figure.figsize': [8, 6], # instead of 4.5, 4.5
	'axes.linewidth': 2,

	'xtick.major.size': 19,
	'ytick.major.size': 19,
	'xtick.minor.size': 10,
	'ytick.minor.size': 10,

	'xtick.major.width': 2,
	'ytick.major.width': 2,
	'xtick.minor.width': 2,
	'ytick.minor.width': 2,

	'xtick.major.pad': 10,
	'ytick.major.pad': 10,
	#'xtick.minor.pad': 14,
	#'ytick.minor.pad': 14,

	'xtick.direction': 'in',
	'ytick.direction': 'in',
   }
matplotlib.rcParams.update(params)

Extract Ophiuchus data from Planck 353 GHz map

https://arxiv.org/pdf/1502.04123.pdf

1. Input Data

We will use Planck HFI 353 GHz map release 3: https://irsa.ipac.caltech.edu/data/Planck/release_3/all-sky-maps/

All sky maps are in HEALPix format, with Nside=2048 (for LFI 70GHz and HFI) with nested ordering.

COORDSYS= 'GALACTIC'
I, Q, U: The signal is given in units of Kcmb for 30-353 GHz.
II, QQ, UU, IQ, IU, QU: The signal in the unit of $K
{cmb}^2$.
Beam: 4.818 arcmins

Unit Convention

From Kelvin to Jansky/sr
$ 1Jy = 10^{-26} W. m^{-2}. Hz ^{-1}$ ; $1W = kg. m^2.s^{−3} = \dfrac{1J}{s} = \dfrac{N.m}{s} = 1V.1A = ...$
sr = square radian
In the Rayleigh-Jeans limit (hν << kT), then

$ 𝑇_{CMB}= \dfrac{𝑆_{\nu} c^2 }{2k_{b} \nu^2 \theta _s} $ [K].
$S_{\nu} = 2 T_{CMB} k_{b} \theta _s \dfrac{ \nu^2}{c^2} $ [Jy/beam, ( or [Jy] assuming it came from a single beam ).

Where:
$𝑆_\nu$ is the flux density [$W m^{-2} Hz^{-1} sr^{-1}$],
c = 299792458 [m/s]
$\nu$ is the frequency [GHz],
${k_b}$ = 1.380649×10−23 [J⋅K−1] is Boltzmann's constant
$\theta_s$ is the size of the source (or the resolution in this case).

example:
beam_sigma = 50" # arcsec
beam_sigma = beam_sigma * $\pi$ / 60 /60 / 180 # radians
$\theta_s$ = beam_area = $2.\pi$.(beam_sigma)$^2$

As long as the source and wavelength are constant, brightness temperature and flux density can be interchanged.

ref0: https://arxiv.org/pdf/1111.1183.pdf
ref1: https://docs.astropy.org/en/stable/api/astropy.units.brightness_temperature.html
ref2: https://astronomy.stackexchange.com/questions/20704/how-are-people-converting-intensities-in-janskys-to-kelvin

In [3]:
kb = 1.380649*1e-23
c=299792458
nu=353*1e9
lamb = c**2 / nu**2   #m
K2Jy = 2*kb*nu**2/c**2 / 1e-26 # [J.K-1.m-2] = [W.s.K-1.m-2] = [W.Hz-1.m-2.K-1]
K2Jy
Out[3]:
3828435046.856711
In [4]:
planck_353 = hp.read_map("HFI_SkyMap_353_2048_R3.01_full.fits",field=(0,1,2),nest=False) #,3,4,5,6,7,8,9))
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/healpy/fitsfunc.py:352: UserWarning: If you are not specifying the input dtype and using the default np.float64 dtype of read_map(), please consider that it will change in a future version to None as to keep the same dtype of the input file: please explicitly set the dtype if it is important to you.
  "If you are not specifying the input dtype and using the default "
NSIDE = 2048
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/healpy/fitsfunc.py:403: UserWarning: No INDXSCHM keyword in header file : assume IMPLICIT
  warnings.warn("No INDXSCHM keyword in header file : " "assume {}".format(schm))
Ordering converted to RING
Ordering converted to RING
Ordering converted to RING
In [5]:
 hp.nside2resol(2048, arcmin = True)
Out[5]:
1.717743205908703
In [6]:
I_stokes = planck_353[0]*K2Jy   
Q_stokes = planck_353[1]*K2Jy
U_stokes = planck_353[2]*K2Jy      
#hits = planck_353[3]
#II_cov = planck_353[4]*1e12    #uK^2
#IQ_cov = planck_353[5]*1e12
#IU_cov = planck_353[6]*1e12
#QQ_cov = planck_353[7]*1e12
#QU_cov = planck_353[8]*1e12
#UU_cov = planck_353[9]*1e12
In [7]:
plt.figure(figsize=(14,10)) ## maps in RING
hp.mollview(I_stokes,coord=["G"],unit='Jy',title="I Stokes",sub=[1,3,1],norm='hist',cmap=cmap)
hp.graticule()
hp.mollview(Q_stokes,coord=["G"],unit='Jy',title="Q Stokes",sub=[1,3,2],norm='hist',cmap=cmap)
hp.graticule()
hp.mollview(U_stokes,coord=["G"],unit='Jy',title="U Stokes",sub=[1,3,3],norm='hist',cmap=cmap)
hp.graticule()
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.

2. Sky regions definition

In [8]:
nside=2048

Ophiuchus_distance=140 #pc
l_ophi=354.0   # Galactic Longitude
b_ophi=17      # Galactic Latitude
delta_l_ophi,delta_b_ophi = 13.,13.

Lupus_distance=140. #pc
l_lupus= 340.0
b_lupus=12.7
delta_l_lupus,delta_b_lupus=12.,12.

Orion_distance = 450 #pc
l_orion = 212.0
b_orion = -16.0
delta_l_orion,delta_b_orion=16.,16.

l_taurus=172.5
b_taurus=-14.5
delta_l_taurus,delta_b_taurus=15.,15.

rot_ophi = [l_ophi,b_ophi]
rot_lupus = [l_lupus,b_lupus]
rot_orion = [l_orion,b_orion]
In [9]:
# uncomment these if you want to see locations
#hp.mollzoom(I_stokes,coord=["G"],title="I Stokes",rot=rot_ophi,norm='hist',cmap=cmap)
#hp.graticule()
#hp.mollzoom(I_stokes,coord=["G"],title="I Stokes",rot=rot_lupus,norm='hist',cmap=cmap)
#hp.graticule()
#hp.mollzoom(I_stokes,coord=["G"],title="I Stokes",rot=rot_orion,norm='hist',cmap=cmap)
#hp.graticule()
In [10]:
plt.figure(figsize=(14,10))
hp.gnomview(I_stokes,coord=["G"],title="I Stokes Ophiuchus",unit='Jy',rot=rot_ophi,sub=[1,3,1],xsize=4000,cmap=cmap,norm='hist')
hp.graticule()
hp.gnomview(I_stokes,coord=["G"],title="I Stokes Lupus",unit='Jy',rot=rot_lupus,sub=[1,3,2],xsize=4000,cmap=cmap,norm='hist')
hp.graticule()
hp.gnomview(I_stokes,coord=["G"],title="I Stokes Orion",unit='Jy',rot=rot_orion,sub=[1,3,3],xsize=2000,cmap=cmap,norm='hist')
hp.graticule()
39.02432778087851 140.9756722191215 -50.97567221912149 50.97567221912149
The interval between parallels is 10 deg 0.00'.
The interval between meridians is 6 deg 0.00'.
39.02432778087851 140.9756722191215 -50.97567221912149 50.97567221912149
The interval between parallels is 10 deg 0.00'.
The interval between meridians is 6 deg 0.00'.
39.02432778087851 140.9756722191215 -50.97567221912149 50.97567221912149
The interval between parallels is 10 deg 0.00'.
The interval between meridians is 6 deg 0.00'.
39.02432778087851 140.9756722191215 -50.97567221912149 50.97567221912149
The interval between parallels is 10 deg 0.00'.
The interval between meridians is 6 deg 0.00'.
39.02432778087851 140.9756722191215 -50.97567221912149 50.97567221912149
The interval between parallels is 10 deg 0.00'.
The interval between meridians is 6 deg 0.00'.
58.335435317857964 121.66456468214203 -31.664564682142032 31.664564682142032
The interval between parallels is 5 deg 0.00'.
The interval between meridians is 4 deg 0.00'.

3. Making masks

In [11]:
def rectangle_mask_region(name,l,b,delta_l,delta_b,nside,write_mask=True,show_mask=True):
    import healpy as hp
    import numpy as np
    # l : Galactic Longitude
    # b : Galactic Latitude

    theta = 90. - b     # Spherical Latitude
    phi = l             # Spherical co-longitude
    
    #calculate rectangle mask
    theta_l,theta_u = theta - delta_b/2., theta + delta_b/2.  # lower, upper 
    phi_l,phi_r = phi - delta_l/2., phi + delta_l/2.          # left, right
    
    v_ll = hp.ang2vec(np.radians(theta_l),np.radians(phi_l))
    v_ul = hp.ang2vec(np.radians(theta_u),np.radians(phi_l))
    v_ur = hp.ang2vec(np.radians(theta_u),np.radians(phi_r))
    v_lr = hp.ang2vec(np.radians(theta_l),np.radians(phi_r))

    vertices = np.vstack((v_ul, v_ll, v_lr, v_ur))
    print('vertices',vertices)
    indexes = hp.query_polygon(nside, vertices)
    print('indexes',indexes)
    
    t, p = np.degrees(hp.pix2ang(nside, indexes))

    NPIX = hp.nside2npix(nside)   ; print(NPIX)
    m = np.arange(NPIX)*0.         #m = np.zeros(hp.nside2npix(nside))
    m[indexes]=1.
    mask = m
    
    if write_mask: hp.write_map('Mask'+name+str(nside)+'_.fits',mask,overwrite=True)
    
    if show_mask:
        hp.mollview(mask,coord=["G"], title="Mollview image RING")
        hp.graticule()
    
    return mask,indexes
In [12]:
mask_ophi,id_ophi=rectangle_mask_region('ophiuchus',l_ophi,b_ophi,delta_l_ophi,delta_b_ophi,nside,write_mask=True,show_mask=False)
mask_lupus,id_lupus=rectangle_mask_region('lupus',l_lupus,b_lupus,delta_l_lupus,delta_b_lupus,nside,write_mask=True,show_mask=False)
mask_orion,id_orion=rectangle_mask_region('orion',l_orion,b_orion,delta_l_orion,delta_b_orion,nside,write_mask=True,show_mask=False)
vertices [[ 0.95994784 -0.21281531  0.18223553]
 [ 0.89532209 -0.19848813  0.39874907]
 [ 0.91702516  0.00800276  0.39874907]
 [ 0.98321747  0.00858041  0.18223553]]
indexes [15085387 15085388 15085389 ... 20582116 20582117 20582118]
50331648
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/healpy/fitsfunc.py:187: FutureWarning: The default dtype of write_map() will change in a future version: explicitly set the dtype if it is important to you
  category=FutureWarning)
vertices [[ 0.89265587 -0.43537736  0.11667074]
 [ 0.85134696 -0.41522966  0.32061299]
 [ 0.91907408 -0.22915091  0.32061299]
 [ 0.96366924 -0.24026973  0.11667074]]
indexes [17067522 17067523 17067524 ... 22228671 22228672 22228673]
50331648
vertices [[-0.8345653  -0.37157241 -0.40673664]
 [-0.9046549  -0.40277831 -0.1391731 ]
 [-0.75858935 -0.63653204 -0.1391731 ]
 [-0.69981642 -0.5872157  -0.40673664]]
indexes [28672546 28672547 28672548 ... 35480326 35480327 35480328]
50331648
In [13]:
hp.gnomview(mask_lupus,coord=["G"],title="mask",rot=rot_lupus,sub=[1,3,3],xsize=600,norm='hist',cmap=cmap)
hp.graticule()

hp.mollview(mask_ophi+mask_lupus+mask_orion,coord=["G"],norm='hist',cmap=cmap)
hp.graticule()
79.52922013126312 100.47077986873688 -10.47077986873688 10.47077986873688
The interval between parallels is 2 deg 30.00'.
The interval between meridians is 1 deg 0.00'.
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.

The standard coordinates are the colatitude 𝜃, 0 at the North Pole, 𝜋/2 at the equator and 𝜋 at the South Pole and the longitude 𝜑 between 0 and 2𝜋 eastward, in a Mollview projection, 𝜑 = 0 is at the center and increases eastward toward the left of the map.

$\theta = \pi/2 - dec $
$\phi = ra $

$\theta = 90 - b$
$\phi = l$

4. Calculate selected regions

In [14]:
I_ophiuchus = I_stokes*mask_ophi
Q_ophiuchus = Q_stokes*mask_ophi
U_ophiuchus = U_stokes*mask_ophi

I_orion = I_stokes*mask_orion
Q_orion = Q_stokes*mask_orion
U_orion = U_stokes*mask_orion
In [15]:
plt.figure(figsize=(14,10))
hp.gnomview(I_ophiuchus,title="I Ophiuchus",unit='Jy',rot=rot_ophi,sub=[1,3,1],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
hp.gnomview(Q_ophiuchus,title="Q Ophiuchus",unit='Jy',rot=rot_ophi,sub=[1,3,2],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
hp.gnomview(U_ophiuchus,title="U Ophiuchus",unit='Jy',rot=rot_ophi,sub=[1,3,3],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
81.5932482582206 98.40675174177939 -8.406751741779399 8.406751741779399
The interval between parallels is 1 deg 0.00'.
The interval between meridians is 1 deg 0.00'.
81.5932482582206 98.40675174177939 -8.406751741779399 8.406751741779399
The interval between parallels is 1 deg 0.00'.
The interval between meridians is 1 deg 0.00'.
81.5932482582206 98.40675174177939 -8.406751741779399 8.406751741779399
The interval between parallels is 1 deg 0.00'.
The interval between meridians is 1 deg 0.00'.
81.5932482582206 98.40675174177939 -8.406751741779399 8.406751741779399
The interval between parallels is 1 deg 0.00'.
The interval between meridians is 1 deg 0.00'.
81.5932482582206 98.40675174177939 -8.406751741779399 8.406751741779399
The interval between parallels is 1 deg 0.00'.
The interval between meridians is 1 deg 0.00'.
81.5932482582206 98.40675174177939 -8.406751741779399 8.406751741779399
The interval between parallels is 1 deg 0.00'.
The interval between meridians is 1 deg 0.00'.
In [16]:
#hp.mollzoom(I_ophiuchus,title="I Stokes",rot=rot_ophi,cmap=cmap,norm='hist')
#hp.graticule()

Polarization degrees:

$p=\dfrac{\sqrt{U^2 + Q^2} }{I}$

Polarization angle is calculated:
$\psi=\dfrac{1}{2}arctan \dfrac{-U}{Q}$
In [17]:
p = np.sqrt(Q_ophiuchus**2 + U_ophiuchus**2) / I_ophiuchus
psi = 0.5*np.arctan(-U_ophiuchus/Q_ophiuchus)      # confuse the sign: - or +

p_orion = np.sqrt(Q_orion**2 + U_orion**2) / I_orion
psi_orion = 0.5*np.arctan(-U_orion/Q_orion)
psi_orion
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in true_divide
  """Entry point for launching an IPython kernel.
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in true_divide
  
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in true_divide
  after removing the cwd from sys.path.
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in true_divide
  """
Out[17]:
array([nan, nan, nan, ..., nan, nan, nan])
In [18]:
# We can calculate just in pixels of Ophiuchus, it will avoid deviding by 0 warning!
#data = np.where(I_ophiuchus!=0)
#p = np.sqrt(Q_ophiuchus[id_ophi]**2 + U_ophiuchus[id_ophi]**2) / I_ophiuchus[id_ophi]
#psi = np.arctan(U_ophiuchus[id_ophi]/Q_ophiuchus[id_ophi])
In [19]:
plt.figure(figsize=(14,10))
hp.gnomview(p,coord=["G"],title="Polarization degrees",rot=rot_ophi,sub=[1,2,1],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
hp.gnomview(psi,coord=["G"],title="Polarization angle",rot=rot_ophi,sub=[1,2,2],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
81.5932482582206 98.40675174177939 -8.406751741779399 8.406751741779399
The interval between parallels is 1 deg 0.00'.
The interval between meridians is 1 deg 0.00'.
81.5932482582206 98.40675174177939 -8.406751741779399 8.406751741779399
The interval between parallels is 1 deg 0.00'.
The interval between meridians is 1 deg 0.00'.
81.5932482582206 98.40675174177939 -8.406751741779399 8.406751741779399
The interval between parallels is 1 deg 0.00'.
The interval between meridians is 1 deg 0.00'.
In [20]:
#m = np.arange(hp.nside2npix(nside))
#plt.plot(m, p)
#plt.xlabel('pixels')
#plt.ylabel('Polarization degrees')

fig, ax = plt.subplots(figsize=(10,7))
ax.hist(p,bins=40,histtype='step',linewidth=4.,color='b',range=[0,0.4],label='Ophiuchus')
ax.hist(p_orion,bins=40,histtype='step',linewidth=4.,color='g',range=[0,0.4],label='Orion')
ax.axvline(x=np.median(p[id_ophi]),linestyle='--',linewidth=4, color='b',label='Median:'+str(round(np.median(p[id_ophi]),3)))
ax.axvline(x=np.mean(p_orion[id_orion]),linestyle='-.',linewidth=4, color='g',label='Median:'+str(round(np.median(p_orion[id_orion]),3)))

plt.xlabel('Polarization degrees')
plt.ylabel('Occurance')
plt.legend()
plt.grid(linestyle='dotted')


fig, ax = plt.subplots(figsize=(10,7))
ax.hist(psi,bins=40,histtype='step',linewidth=4.,color='b',label='Ophiuchus')#,range=[0,0.6])
ax.hist(psi_orion,bins=40,histtype='step',linewidth=4.,color='g',label='Orion')#,range=[0,0.6])
ax.axvline(x=np.median(psi[id_ophi]),linestyle='--',linewidth=4, color='b',label='Median:'+str(round(np.median(psi[id_ophi]),3)))
ax.axvline(x=np.mean(psi_orion[id_orion]),linestyle='--',linewidth=4, color='g',label='Median:'+str(round(np.median(psi_orion[id_orion]),3)))

plt.xlabel('Polarization angle [rad]')
plt.ylabel('Occurance')
plt.legend(loc='lower center')
plt.grid(linestyle='dotted')
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/lib/histograms.py:839: RuntimeWarning: invalid value encountered in greater_equal
  keep = (tmp_a >= first_edge)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/lib/histograms.py:840: RuntimeWarning: invalid value encountered in less_equal
  keep &= (tmp_a <= last_edge)
ophi_vec = hp.pix2vec(nside,id_ophi) # a test in vector

5. 2D real space data maps

These are mollzoom plots to test the location we are going to take the data, they are just different in normalization color.

hp.mollzoom(I_stokes,coord=["G"],title="I Ophiuchus",unit='Jy',rot=rot_ophi,cmap=cmap,norm='hist')
hp.mollzoom(I_stokes,coord=["G"],title="I Ophiuchus",unit='Jy',rot=rot_ophi,cmap=cmap)

In [21]:
data = I_ophiuchus[id_ophi]   # a linear dimension data
np.size(data)
np.size(id_ophi)
Out[21]:
197153
In [22]:
histogram = plt.hist(data,bins=1000,range=[0,1e8])
In [23]:
NSIDE = 2048
hp.nside2resol(NSIDE, arcmin = True)
Out[23]:
1.717743205908703
In [24]:
N=np.size(data) #=197153 this is the number of pixels in a linear dimension
                # since we are using lots of FFTs this should be a factor of 2^N
pix_size = hp.nside2resol(NSIDE, arcmin = True) #1.7177 # size of a pixel in arcminutes
planck_beam=4.8 #arc

## variables to set up the map plots; we have 197153 pixels, each pixel has 1.7177 arcminutes
X_width = N*pix_size/60. # horizontal map width in degrees , /60 = convert arcminutes to degrees
Y_width = N*pix_size/60. # vertical map width in degrees
                                #(we can devided by 2 for fast now) 
print('X_width',X_width,'Y_width',Y_width)
X_width 5644.3037712419755 Y_width 5644.3037712419755
In [25]:
import numpy as np
import matplotlib.pylab as plt
from pylab import arange, show, cm
cmap = cm.bwr
cmap.set_under('w')

from astropy.io import fits
hdulist = fits.open('HFI_SkyMap_353_2048_R3.01_full.fits')  
hdulist.info()
#hdulist[1].header
Filename: HFI_SkyMap_353_2048_R3.01_full.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       4   ()      
  1  FREQ-MAP      1 BinTableHDU     75   50331648R x 10C   [E, E, E, J, E, E, E, E, E, E]   
# Maps in NESTED plt.figure(figsize=(14,10)) hp.mollview(hdulist[1].data['I_STOKES']*K2Jy,coord=["G"],norm='hist',cmap=cmap,sub=[1,3,1]) # in Galactic coordinate, no NEST or RING? hp.mollview(hdulist[1].data['Q_STOKES']*K2Jy,norm='hist',cmap=cmap,sub=[1,3,2]) hp.mollview(hdulist[1].data['U_STOKES']*K2Jy,norm='hist',cmap=cmap,sub=[1,3,3])
In [26]:
## Test Interpolate values from healpy maps
#I_stokes = I_stokes*K2Jy  # These maps in RING format, so if we use this we have to change hpx = HEALPix(nside=nside, order='RING', frame=Galactic())
#Q_stokes = I_stokes*K2Jy  #
#U_stokes = I_stokes*K2Jy  #

## Interpolate values from FITS file
I_stokes = hdulist[1].data['I_STOKES']*K2Jy   # These data maps in NESTED format
Q_stokes = hdulist[1].data['Q_STOKES']*K2Jy   #
U_stokes = hdulist[1].data['U_STOKES']*K2Jy   #
## we will use astropy interpolated bilinear method; another method as color excess method: https://core.ac.uk/reader/4881074
In [27]:
# Set up the HEALPix projection of the data: the 353GHz map is in the Galactic coordinate
from astropy_healpix import HEALPix
from astropy.coordinates import Galactic
nside = hdulist[1].header['NSIDE']      # 2048
order = hdulist[1].header['ORDERING']   # NESTED
hpx = HEALPix(nside=nside, order=order, frame=Galactic())   # Make a map in Galactic frame ; NEST

5.1. In Galactic coordinate (longitude l,lattitude b)

In [28]:
# Sample a X_width x Y_width grid which is coressponding to l,b, it depends on the data region.
from astropy import units as u
l=np.linspace( l_ophi - delta_l_ophi/2., l_ophi + delta_l_ophi/2., int(X_width)) * u.deg
b=np.linspace(b_ophi - delta_b_ophi/2., b_ophi + delta_l_ophi /2., int(Y_width)) * u.deg
l_grid,b_grid=np.meshgrid(l,b)

# Set up Astropy coordinate objects, the coordinate we want to see the objects
from astropy.coordinates import SkyCoord
coords = SkyCoord(l_grid.ravel(), b_grid.ravel(), frame='galactic')
In [29]:
Imap = hpx.interpolate_bilinear_skycoord(coords, I_stokes)
Qmap = hpx.interpolate_bilinear_skycoord(coords, Q_stokes)
Umap = hpx.interpolate_bilinear_skycoord(coords, U_stokes)

Imap_gal = np.reshape(Imap, ( int(X_width), int(Y_width) ) )[::-1,::-1]  #swap x,y axis
Qmap_gal = np.reshape(Qmap, ( int(X_width), int(Y_width) ) )[::-1,::-1]
Umap_gal = np.reshape(Umap, ( int(X_width), int(Y_width) ) )[::-1,::-1]
In [30]:
#histogram = plt.hist(I_map,bins=1000,range=[0,1e8])
In [31]:
from matplotlib.colors import LogNorm
from mpl_toolkits.axes_grid1 import make_axes_locatable

ex=[l.value.min(), l.value.max(), b.value.min(), b.value.max() ]
plt.figure(figsize=(10,8) )
im = plt.imshow(Imap_gal, extent=ex , norm=LogNorm(), cmap=cmap) #,origin='left' )
plt.colorbar(im,pad=0.05)
plt.xlabel('Galactic longtitude $l$ [deg]')
plt.ylabel('Galactic latitude $b$ [deg]')
plt.show()
In [32]:
from mpl_toolkits.axes_grid1 import make_axes_locatable

fig = plt.figure(figsize=(16, 12))
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=None)

ax1 = fig.add_subplot(121)
plt.xlabel(r"Galactic longtitude $l$ [deg]",fontsize="20")
plt.ylabel(r"Galactic latitude $b$ [deg]",fontsize="20")
im1 = ax1.imshow(Qmap_gal, extent=ex , norm=LogNorm(), cmap=cmap) #,origin='left' )
divider = make_axes_locatable(ax1)
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(im1, cax=cax, orientation='vertical')


ax2 = fig.add_subplot(122)
plt.xlabel(r"Galactic longtitude $l$ [deg]",fontsize="20")
plt.ylabel(r"Galactic latitude $b$ [deg]",fontsize="20")
im2 = ax2.imshow(Umap_gal,extent=ex, norm=LogNorm(), cmap=cmap)
divider = make_axes_locatable(ax2)
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(im2, cax=cax, orientation='vertical')
Out[32]:
<matplotlib.colorbar.Colorbar at 0x16248fc88>
In [33]:
outfile_I='I_stokes_Ophichus_Galactic_coord.fits'   # HDU = Header Data Unit
outfile_Q='Q_stokes_Ophichus_Galactic_coord.fits'   # HDU = Header Data Unit
outfile_U='U_stokes_Ophichus_Galactic_coord.fits'   # HDU = Header Data Unit

hdu=fits.PrimaryHDU(Imap_gal)
hdu.writeto(outfile_I, clobber=True)

hdu=fits.PrimaryHDU(Qmap_gal)
hdu.writeto(outfile_Q, clobber=True)

hdu=fits.PrimaryHDU(Umap_gal)
hdu.writeto(outfile_U, clobber=True)
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 2.0 and will be removed in a future version. Use argument "overwrite" instead. [__main__]
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 2.0 and will be removed in a future version. Use argument "overwrite" instead. [__main__]
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 2.0 and will be removed in a future version. Use argument "overwrite" instead. [__main__]

5.2. In equatorial coordinate (RA,DEC)

In [34]:
from astropy import units as u
from astropy.coordinates import SkyCoord
c = SkyCoord(l_ophi, b_ophi, unit='deg', frame='galactic')
c.icrs
Out[34]:
<SkyCoord (ICRS): (ra, dec) in deg
    (247.1758789, -23.6727017)>
In [35]:
ra_ophi = c.icrs.ra.hms 
ra_ophi
Out[35]:
hms_tuple(h=16.0, m=28.0, s=42.210936665427994)
In [36]:
# Sample a X_width x Y_width grid which is coressponding to RA,DEC of region.
ra_ophi = c.icrs.ra.degree
dec_ophi = c.icrs.dec.degree
ra = np.linspace(ra_ophi - delta_l_ophi/2., ra_ophi + delta_l_ophi /2., int(X_width)) * u.deg
dec = np.linspace(dec_ophi - delta_b_ophi /2., dec_ophi + delta_l_ophi /2., int(Y_width)) * u.deg
ra_grid, dec_grid = np.meshgrid(ra, dec)

# Set up Astropy coordinate objects
from astropy.coordinates import SkyCoord
coords = SkyCoord(ra_grid.ravel(), dec_grid.ravel(), frame='icrs')

# Interpolate values
I_map = hpx.interpolate_bilinear_skycoord(coords, I_stokes)
Q_map = hpx.interpolate_bilinear_skycoord(coords, Q_stokes)
U_map = hpx.interpolate_bilinear_skycoord(coords, U_stokes)

Imap_equa = np.reshape(I_map, ( int(X_width), int(Y_width) ) )[::-1,::-1]   #.swapaxes(-1,-1)#[...,::-1] #(-1,0)[...,::-2] #.swapaxes(-1,0)
Qmap_equa = np.reshape(Q_map, ( int(X_width), int(Y_width) ) )[::-1,::-1]
Umap_equa = np.reshape(U_map, ( int(X_width), int(Y_width) ) )[::-1,::-1]
In [37]:
from matplotlib.colors import LogNorm
ex=[ ra.value.max(),ra.value.min(),dec.value.min(), dec.value.max()]
plt.figure(figsize=(10,8))
im = plt.imshow(Imap_equa, extent=ex,norm=LogNorm(),cmap=cmap)#,origin='right')
plt.colorbar(im)

plt.xlabel('Right ascension [deg]')
plt.ylabel('Declination [deg]')
plt.show()
In [38]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
ex=[ ra.value.max(),ra.value.min(),dec.value.min(), dec.value.max()]
fig = plt.figure(figsize=(16, 12))
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=None)

ax1 = fig.add_subplot(121)
plt.xlabel('Right ascension [deg]')
plt.ylabel('Declination [deg]')
im1 = ax1.imshow(Qmap_equa, extent=ex , norm=LogNorm(), cmap=cmap) #,origin='left' )
divider = make_axes_locatable(ax1)
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(im1, cax=cax, orientation='vertical')


ax2 = fig.add_subplot(122)
plt.xlabel('Right ascension [deg]')
plt.ylabel('Declination [deg]')
im2 = ax2.imshow(Umap_equa,extent=ex, norm=LogNorm(), cmap=cmap)
divider = make_axes_locatable(ax2)
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(im2, cax=cax, orientation='vertical')
Out[38]:
<matplotlib.colorbar.Colorbar at 0x162d8bef0>
In [39]:
outfile_I='I_stokes_Ophichus_RADEC_coord.fits'   # HDU = Header Data Unit
outfile_Q='Q_stokes_Ophichus_RADEC_coord.fits'   # HDU = Header Data Unit
outfile_U='U_stokes_Ophichus_RADEC_coord.fits'   # HDU = Header Data Unit


hdu=fits.PrimaryHDU(Imap_equa)
hdu.writeto(outfile_I, clobber=True)

hdu=fits.PrimaryHDU(Qmap_equa)
hdu.writeto(outfile_Q, clobber=True)

hdu=fits.PrimaryHDU(Umap_equa)
hdu.writeto(outfile_U, clobber=True)
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 2.0 and will be removed in a future version. Use argument "overwrite" instead. [__main__]
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 2.0 and will be removed in a future version. Use argument "overwrite" instead. [__main__]
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 2.0 and will be removed in a future version. Use argument "overwrite" instead. [__main__]

6. Test written fits data

In [40]:
infile='I_stokes_Ophichus_Galactic_coord.fits'
In [41]:
hdu_list=fits.open(infile)
hdu_list.info()
Filename: I_stokes_Ophichus_Galactic_coord.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       6   (5644, 5644)   float64   
In [42]:
I_ophi_data = hdu_list[0].data
print(type(I_ophi_data))
print(I_ophi_data.shape)
<class 'numpy.ndarray'>
(5644, 5644)
In [43]:
plt.imshow(I_ophi_data,norm=LogNorm(),cmap=cmap)
Out[43]:
<matplotlib.image.AxesImage at 0x1662e8d68>
In [ ]:
 

7. References

[0] https://arxiv.org/pdf/1502.04123.pdf

[1] UNDERSTANDING POLARIZED DUST EMISSION FROM ρ OPHIUCHI A IN LIGHT OF GRAIN ALIGNMENT AND DISRUPTION BY RADIATIVE TORQUES https://arxiv.org/pdf/2007.10621.pdf

[2] Planck Early Results. XXV. Thermal dust in nearby molecular clouds https://arxiv.org/pdf/1101.2037.pdf

[3] https://arxiv.org/pdf/1909.04862.pdf

Appendix

## LIC method: update in healpy 1.40 import numpy as np import matplotlib.pyplot as plt from matplotlib.cm import get_cmap from matplotlib.colors import ListedColormap import healpy as hp cmap_colors = get_cmap('binary', 256)(np.linspace(0, 1, 256)) #cmap_colors[..., 3] = 0.4  # Make colormap partially transparent cmap = cmap # ListedColormap(cmap_colors) m = hp.read_map('HFI_SkyMap_353_2048_R3.01_full.fits', (0, 1, 2), verbose=False) I, Q, U = hp.smoothing(m, np.deg2rad(5)) lic = hp.line_integral_convolution(Q, U) lic = hp.smoothing(lic, np.deg2rad(0.5)) hp.mollview(np.log(1 + np.sqrt(Q**2 + U**2) * 100), cmap='inferno', cbar=False) hp.mollview(lic, cmap=cmap, cbar=False, reuse_axes=True, title='WMAP K') #plt.savefig('wmapk.png', dpi=150)
In [44]:
# Set up Astropy coordinate objects
#from astropy.coordinates import SkyCoord
#coord = SkyCoord( frame='galactic', l="354.", b="17.",unit='deg')
In [45]:
from math import *
import numpy as np

def ga2equ(ga):
    """
    Convert Galactic to Equatorial coordinates (J2000.0)
    (use at own risk)
    
    Input: [l,b] in decimal degrees
    Returns: [ra,dec] in decimal degrees
    
    Source: 
    - Book: "Practical astronomy with your calculator" (Peter Duffett-Smith)
    - Wikipedia "Galactic coordinates"
    
    Tests (examples given on the Wikipedia page):
    >>> ga2equ([0.0, 0.0]).round(3)
    array([ 266.405,  -28.936])
    >>> ga2equ([359.9443056, -0.0461944444]).round(3)
    array([ 266.417,  -29.008])
    """
    l = radians(ga[0])
    b = radians(ga[1])

    # North galactic pole (J2000) -- according to Wikipedia
    pole_ra = radians(192.859508)
    pole_dec = radians(27.128336)
    posangle = radians(122.932-90.0)
    
    # North galactic pole (B1950)
    #pole_ra = radians(192.25)
    #pole_dec = radians(27.4)
    #posangle = radians(123.0-90.0)
    
    ra = atan2( (cos(b)*cos(l-posangle)), (sin(b)*cos(pole_dec) - cos(b)*sin(pole_dec)*sin(l-posangle)) ) + pole_ra
    dec = asin( cos(b)*cos(pole_dec)*sin(l-posangle) + sin(b)*sin(pole_dec) )
    
    return np.array([degrees(ra), degrees(dec)])

ra,dec=ga2equ([354.,17.])
ra,dec
Out[45]:
(247.17582815053345, -23.67270697245831)